#### Load Packages
library(languageserver)
library(readxl)
library(sandwich)
library(lmtest)
library(tidyverse)
library(haven)
library(stargazer)
library(foreign)
library(manifestoR)
library(countrycode)
library(lubridate)
library(lfe)
library(xtable)
library(sjPlot)
library(sjmisc)
library(ggeffects)
library(coefplot)
library(broom)
library(devtools)



##
rm(list=ls())
options(scipen = 999)

getwd()

################################################################################
# Load Data

full_dat <- readRDS('data/WEP_2024_data.RDS') 

# define data without missing positions (restricted sample)
na_dat <- full_dat %>%  drop_na(CAP_immig_kimfor_na, EP_immig_kimfor_na) 

####################################################
# Create Table 1:

ols_m1_change <- felm(kimfor_change ~    radical_parliament  + factor(mwc) +
              no_cabparties + unempl_cab_begin 
              | iso2c + year |0|iso2c,data=full_dat)

ols_m3_change <- felm(kimfor_change ~    radical_parliament  + factor(mwc) +
              no_cabparties + unempl_cab_begin
               | iso2c + year |0|iso2c,data=na_dat)

ols_m2_change <- felm(log_change ~    radical_parliament  + factor(mwc) +
              no_cabparties + unempl_cab_begin 
               |  iso2c + year |0|iso2c ,data=full_dat)

ols_m4_change <- felm(log_change ~    radical_parliament  + factor(mwc) +
              no_cabparties + unempl_cab_begin 
                | iso2c + year |0|iso2c ,data=na_dat)

stargazer( ols_m1_change, ols_m3_change, ols_m2_change, ols_m4_change,
           title="Table 1",
           digits = 3,
           ci.level=0.95, 
           add.lines = list(c("Fixed Effects", "Country+Year", "Country+Year", "Country+Year", "Country+Year") ),
           omit.stat = c("rsq", "f","ser"),
           omit = "as.factor",
           type = 'text',
           no.space=TRUE
           , out= c("output/table_1.tex",  "output/table_1.html")
            )


####################################################
# Create Table 2:
ols_m1_change <- felm(kimfor_change ~   radical_parliament + radical_votes_total +  factor(mwc) +
              no_cabparties + unempl_cab_begin 
                 | iso2c + year |0|iso2c,data=full_dat)

ols_m2_change <- felm(log_change ~    radical_parliament + radical_votes_total +  factor(mwc) +
              no_cabparties + unempl_cab_begin 
               |  iso2c + year |0|iso2c ,data=full_dat)

ols_m3_change <- felm(kimfor_change ~    radical_parliament + radical_votes_total +  factor(mwc) +
              no_cabparties + unempl_cab_begin 
                      | iso2c + year |0|iso2c,data=na_dat)

ols_m4_change <- felm(log_change ~   radical_parliament + radical_votes_total  +  factor(mwc) +
              no_cabparties + unempl_cab_begin 
                      | iso2c + year |0|iso2c ,data=na_dat)

stargazer( ols_m1_change, ols_m3_change, ols_m2_change, ols_m4_change,
           title="Table 2",
           digits = 3,
           ci.level=0.95, 
           add.lines = list(c("Fixed Effects", "Country+Year", "Country+Year", "Country+Year", "Country+Year") ),
           omit.stat = c("rsq", "f","ser"),
           omit = "as.factor",
           type = 'text',
           no.space=TRUE
           , out= c("output/table_2.tex","output/table_2.html")                                  
            )



########################
# Create Table 3

# Interaction with Performance
ols_m1_perf <- felm(kimfor_change ~ radical_parliament * lost_votes + factor(mwc) +
                        no_cabparties + unempl_cab_begin
                      | iso2c + year | 0 | iso2c, data = full_dat)

ols_m3_perf <- felm(kimfor_change ~ radical_parliament * lost_votes + factor(mwc) +
                        no_cabparties + unempl_cab_begin
                      | iso2c + year | 0 | iso2c, data = na_dat)

# Interaction with Left Cabinet
ols_m1_rile <- felm(kimfor_change ~ radical_parliament * left_cab + factor(mwc) +
                        no_cabparties + unempl_cab_begin
                      | iso2c + year | 0 | iso2c, data = full_dat)

ols_m3_rile <- felm(kimfor_change ~ radical_parliament * left_cab + factor(mwc) +
                        no_cabparties + unempl_cab_begin
                      | iso2c + year | 0 | iso2c, data = na_dat)

# Interaction with Long Cabinet Bargaining Duration
ols_m1_bargain <- felm(kimfor_change ~ radical_parliament * long_cab + factor(mwc) +
                           no_cabparties + unempl_cab_begin
                      | iso2c + year | 0 | iso2c, data = full_dat)

ols_m3_bargain <- felm(kimfor_change ~ radical_parliament * long_cab + factor(mwc) +
                           no_cabparties + unempl_cab_begin
                      | iso2c + year | 0 | iso2c, data = na_dat)

# Interaction with Minority Government
ols_m1_minority <- felm(kimfor_change ~ radical_parliament * minority + factor(mwc) +
                            no_cabparties + unempl_cab_begin
                      | iso2c + year | 0 | iso2c, data = full_dat)

ols_m3_minority <- felm(kimfor_change ~ radical_parliament * minority + factor(mwc) +
                            no_cabparties + unempl_cab_begin
                      | iso2c + year | 0 | iso2c, data = na_dat)

# Combined table for all interactions
stargazer(ols_m1_perf, ols_m3_perf, ols_m1_rile, ols_m3_rile,
          ols_m1_bargain, ols_m3_bargain, ols_m1_minority, ols_m3_minority,
          title = "Table 3", digits = 2,
          dep.var.labels = c("Change in Coalition Agreement Position"),
          add.lines = list(c("Fixed Effects", "Country+Year", "Country+Year", "Country+Year", "Country+Year",
                             "Country+Year", "Country+Year", "Country+Year", "Country+Year")),
          omit.stat = c("rsq", "f", "ser"), omit = "as.factor", type = 'text',
          no.space = TRUE, out = c("output/table_3.tex", "output/table_3.html"))





########################
# Create Table 4

# Define the variables 
variables <- c("immig_ln", "surplus", "minority", "elec_perf_1", "rel_cab_dur_own", "enp", "cab_rile", "east")

test_results <- purrr::map_df(variables, function(var) {
  test <- t.test(reformulate("radical_parliament", response = var), data = full_dat)
  tibble(
    variable = var,
    estimate1 = round(as.numeric(test$estimate[1]), 2),
    estimate2 = round(as.numeric(test$estimate[2]), 2),
    p_value = round(as.numeric(test$p.value), 2)
  )
})

test_results <- test_results %>% 
  mutate(
    balance = ifelse(p_value < 0.05, "No", "Yes")  # balance column based on p_value
      ) %>%
  select(-p_value)

# Display the table
print(test_results)

print(xtable(test_results), type = "latex", file = "output/table_4.tex")


########################
# Create Figure 2
# For Figure 2a: Visualize trends for 'kimfor_full'
trends_kimfor_full <- full_dat %>% 
  select(cab_id_erdda, radical_parliament, contains('kimfor_full')) %>%
  pivot_longer(cols = c("EP_immig_kimfor_full_lag", "EP_immig_kimfor_full", "CAP_immig_kimfor_full"), 
               names_to = "pos_type", 
               values_to = "pos_immig") %>% 
  mutate(
    pos_type = factor(
      pos_type, 
      levels = c("EP_immig_kimfor_full_lag", "EP_immig_kimfor_full", 'CAP_immig_kimfor_full')
    ),
    election = case_when(
      pos_type == "EP_immig_kimfor_full_lag"  ~ -1,
      pos_type == "EP_immig_kimfor_full"  ~ 0,
      pos_type == "CAP_immig_kimfor_full"  ~ 1
    )
  )

p1 <- ggplot(trends_kimfor_full, aes(x=election, y=pos_immig, color=as.factor(radical_parliament))) +
  geom_point(shape=1) +
  geom_smooth() +
  theme_bw() +
  theme(legend.position="bottom") +
  ylab("Immigration Position (Kim/Fording)") +
  scale_x_continuous(
    limits = c(-1, 1), breaks = c(-1, 0, 1), labels=c("EP (t-1)", "EP (t0)", "CAP")
  ) +
  labs(colour= "RRP in Parliament")

p1
ggsave('output/figure_2a.png', p1, width = 20, height = 15, units = "cm")

# For Figure 2b: Visualize trends for 'log_full'
trends_log_full <- full_dat %>%
  select(cab_id_erdda, radical_parliament, contains('log_full')) %>%
  pivot_longer(cols = c("EP_immig_log_full_lag", "EP_immig_log_full", "CAP_immig_log_full"), 
               names_to = "pos_type", 
               values_to = "pos_immig") %>%
  mutate(
    pos_type = factor(
      pos_type, 
      levels = c("EP_immig_log_full_lag", "EP_immig_log_full", 'CAP_immig_log_full')
    ),
    election = case_when(
      pos_type == "EP_immig_log_full_lag"  ~ -1,
      pos_type == "EP_immig_log_full"  ~ 0,
      pos_type == "CAP_immig_log_full"  ~ 1
    )
  )

p2 <- ggplot(trends_log_full, aes(x=election, y=pos_immig, color=as.factor(radical_parliament))) +
  geom_point(shape=1) +
  geom_smooth() +
  theme_bw() +
  theme(legend.position="bottom") +
  ylab("Immigration Position (Log-Scale)") +
  scale_x_continuous(
    limits = c(-1, 1), breaks = c(-1,0,1), labels=c("EP (t-1)", "EP (t0)", "CAP")
  ) +
  labs(colour= "RRP in Parliament")

p2
ggsave('output/figure_2b.png', p2, width = 20, height = 15, units = "cm")

